perm filename SUBSLM.OLD[DRW,LCS] blob sn#108366 filedate 1974-12-13 generic text, type T, neo UTF8
00100		SUBROUTINE SS
00200		COMMON X(100),Y(100),N,X1(512),Y1(512),S(100),K
00300		COMMON/SLM/A(2,200),B(2,200),C(2,200),I(200),I1(200),
00350		1 I2(200),KT,TT
00360		DATA L1/1/,L2/2/,L0/0/
00400	C	FIND MAX AND MIN POINTS AND SET SLOPES
00500		ICOUNT=0
00600	710	DO 100 K=2,N-1
00610		KP=K-1
00620		KQ=K+1
00700		IF(Y(K).GT.Y(KP).AND.Y(K).GT.Y(KQ)) GO TO 20
00800		IF(Y(K).LT.Y(KP).AND.Y(K).LT.Y(KQ)) GO TO 20
00900		IF(X(K).GT.X(KP).AND.X(K).GT.X(KQ)) GO TO 30
01000		IF(X(K).LT.X(KP).AND.X(K).LT.X(KQ)) GO TO 30
01100		GO TO 100
01200	20	S(K)=0
01300		ICOUNT=ICOUNT+1
01400		I(ICOUNT)=K
01500		GO TO 100
01600	30	S(K)=1001
01700		ICOUNT=ICOUNT+1
01800		I(ICOUNT)=K
01900	100	CONTINUE
02000		ICOUNT=ICOUNT+1
02100		I(ICOUNT)=N
02200	C	FIND POINTS SUCH THAT X(K)=X(K+1) OR Y(K)=Y(K+1)
02300		DO 900 K=2,N-1
02400		IF(Y(K).EQ.Y(K+1)) I2(K)=1
02500		IF(X(K).EQ.X(K+1)) I2(K)=1
02600	900	CONTINUE
02700	C	SET THE SLOPES BETWEEN
02800		IA=1
02900		IX=0
03000	110	IF(IX.EQ.ICOUNT) GO TO 200
03100		IX=IX+1
03200		IA=IA+1
03300	120	IF(IA.EQ.I(IX)) GO TO 110
03310		KP=IA+1
03400		S(IA)=(Y(KP)-Y(IA-1))/(X(KP)-X(IA-1))
03500		IA=KP
03600		GO TO 120
03700	200	CONTINUE
03800	C	SET BEGIN AND END SLOPES
03900	CC	SZ=1002001.
04000		K=1
04100		CALL S101
04200		CALL S102(X(1),Y(1),X(2),Y(2),S(1),S(2))
04300		IF(S(K).EQ.-1001.) S(K)=1001.
04400		IF(S(K+1).EQ.-1001.) S(K+1)=1001.
04500		K=N-1
04600		CALL S101
04700		CALL S102(X(N),Y(N),X(N-1),Y(N-1),S(N),S(N-1))
04800		IF(S(K).EQ.-1001.) S(K)=1001.
04900		IF(S(K+1).EQ.-1001.) S(K+1)=1001.
05000	C	RESET THE SLOPES FOR STRAIGHT LINES
05100		DO 610 K=1,N-2
05200		U1=SLM2(L0)
05300		U2=SLM2(L1)
05400		IF(ABS(U1-U2).GT..0001) GO TO 610
05500		I1(K)=1
05600		I1(K+1)=1
05700		S(K)=U1
05800		S(K+1)=U1
05900		S(K+2)=U1
06000	610	CONTINUE
06100	C	ADD POINTS
06200		K1=N-1
06300		DO 300 K=1,K1
06400		IF(I1(K).EQ.1) GO TO 300
06500		IF(I2(K).EQ.1) GO TO 300
06600		FLG=0.
06700		IF(S(K).EQ.0..AND.S(K+1).EQ.0.) FLG=1.
06800		IF(S(K).EQ.1001..AND.S(K+1).EQ.1001.) FLG=2.
06900		IF(FLG.EQ.1..OR.FLG.EQ.2.) GO TO 205
07000	CC	SZ=1002001.
07100		CALL S101
07110		KP=K+1
07200		U=(Y(K)-Y(KP)-S(K)*(X(K)-X(KP)))/(S(KP)-S(K))
07300		AX=X(KP)+U
07400		AY=Y(KP)+S(KP)*U
07500		XA=X(K)
07600	1610	XB=X(KP)
07700		IF(XA.LE.XB) GO TO 202
07800		XA=X(KP)
07900		XB=X(K)
08000	202	YA=Y(K)
08100		YB=Y(KP)
08200		IF(YA.LE.YB) GO TO 204
08300		YA=Y(KP)
08400		YB=Y(K)
08500	204	CONTINUE
08600		IF(AX.GE.XA.AND.AX.LE.XB.AND.
08700	     1	AY.GE.YA.AND.AY.LE.YB) GO TO 290
08800	205	K1=K1+1
08900		DO 210 K2=K+1,N
09000		K3=N+K+1-K2
09010		KP=K3+1
09100		X(KP)=X(K3)
09200		Y(KP)=Y(K3)
09300		I1(KP)=I1(K3)
09400		I2(KP)=I2(K3)
09500	210	S(KP)=S(K3)
09600		IF(FLG.EQ.1..OR.FLG.EQ.2.) GO TO 280
09700		N=N+1
09710		KP=K+1
09720		KQ=K+2
09800		Z0=X(K)**2-X(KQ)**2
09900		Z1=Y(K)-Y(KQ)-S(K)*(X(K)-X(KQ))-((X(K)+X(KQ))/2.
10000	     1	-X(K))*(S(K)-S(KQ))
10100	     	Z2=X(K)**3-X(KQ)**3-3*(X(K)-X(KQ))*X(K)**2
10200	     1	-1.5*(X(K)+X(KQ))*Z0+3*X(K)*Z0
10300		AZ=Z1/Z2
10400		BZ=(S(K)-S(KQ)-3*(X(K)**2-X(KQ)**2)*AZ)/
10500	     1	(2*(X(K)-X(KQ)))
10600		CZ=S(K)-3*AZ*X(K)**2-2*BZ*X(K)
10700		DZ=Y(K)-AZ*X(K)**3-BZ*X(K)**2-CZ*X(K)
10800		X(KP)=-BZ/(3.*AZ)
10900		Y(KP)=AZ*X(KP)**3+BZ*X(KP)**2+CZ*X(KP)+DZ
11000		S(KP)=3*AZ*X(KP)**2+2*BZ*X(KP)+CZ
11100		IF(Y(KP).LT.YB.AND.Y(KP).GT.YA) GO TO 290
11200		X(KP)=(X(KQ)+X(K))/2.
11300		Y(KP)=(Y(KQ)+Y(K))/2.
11400		S(KP)=0.
11500		GO TO 290
11600	280	N=N+1
11700		IF(FLG.EQ.2.) GO TO 282
11800		S(K+1)=2.*(Y(K+2)-Y(K))/(X(K+2)-X(K))
11900		GO TO 284
12000	282	S(K+1)=(Y(K+2)-Y(K))/(2.*(X(K+2)-X(K)))
12100	284	X(K+1)=(X(K+2)+X(K))/2.
12200		Y(K+1)=(Y(K+2)+Y(K))/2.
12300	290	IF(S(K).EQ.-1001.) S(K)=1001.
12400		IF(S(K+1).EQ.-1001.) S(K+1)=1001.
12500	300	CONTINUE
12600	C	FIT THE POINTS WITH N-1 CURVES
12700	305	DO 400 K=1,N-1
12750		KP=K+1
12800		IF(I1(K).EQ.0) GO TO 310
12900		A(1,K)=0.
13000		A(2,K)=0.
13100		B(1,K)=X(KP)-X(K)
13200		B(2,K)=Y(KP)-Y(K)
13300		C(1,K)=X(K)
13400		C(2,K)=Y(K)
13500		GO TO 400
13600	310	CALL S101
13700		B(1,K)=2*(Y(KP)-Y(K)-S(KP)*(X(KP)-X(K)))/
13800	     1	(S(K)-S(KP))
13900		B(2,K)=S(K)*B(1,K)
14000		A(1,K)=X(KP)-X(K)-B(1,K)
14100		A(2,K)=Y(KP)-Y(K)-S(K)*B(1,K)
14200		C(1,K)=X(K)
14300		C(2,K)=Y(K)
14400	400	CONTINUE
14500	C	CALCULATE 512 POINTS
14600		M=N-1
14700		R=M/511.
14800		T=-R
14900		DO 500 L=1,511
15000		T=T+R
15100		KT=T+1
15200		KU=T
15300		TT=T-KU
15400		X1(L)=SLM3(L1)
15500	500	Y1(L)=SLM3(L2)
15600		X1(512)=X(N)
15700		Y1(512)=Y(N)
15800	600	RETURN
15900		END
16000	
16100	
16200	C	GIVEN TWO POINTS (XX1,YY1), (XX2,YY2) AND A SLOPE SS2
16300	C	THIS ROUTINE FINDS THE SLOPE SS1
16400		SUBROUTINE S102(XX1,YY1,XX2,YY2,SS1,SS2)
16500		LF=0
16600		IF(SS2.NE.0.) GO TO 5
16700		IF((XX1.LT.XX2.AND.YY1.LT.YY2).OR.
16800	     1	(XX2.LT.XX1.AND.YY2.LT.YY1)) SS1=1001.
16900		IF((XX1.LT.XX2.AND.YY1.GT.YY2).OR.
17000	     1	(XX2.LT.XX1.AND.YY2.GT.YY1)) SS1=-1001.
17100		GO TO 100
17200	5	T=(YY1-YY2)/SS2
17300		X=XX2+T
17400		Y=YY1
17500		XMAX=XX2
17600		XMIN=XX1
17700		IF(XX2.GE.XX1) GO TO 10
17800		XMAX=XX1
17900		XMIN=XX2
18000	10	IF(X.GT.XMIN.AND.X.LT.XMAX) GO TO 20
18100		T=XX1-XX2
18200		Y=YY2+SS2*T
18300		X=XX1
18400		LF=1
18500	20	IF(LF.EQ.1) GO TO 30
18600		D2=XX1-X
18700		D1=X-XX2
18800		R=D2/D1
18900		XX=XX2
19000		YY=(YY2+R*YY1)/(1.+R)
19100		GO TO 40
19200	30	D2=YY1-Y
19300		D1=Y-YY2
19400		R=D2/D1
19500		YY=YY2
19600		XX=(XX2+R*XX1)/(1.+R)
19700	40	SS1=(YY-YY1)/(XX-XX1)
19800	100	RETURN
19900		END
20000	
20100		SUBROUTINE S101
20200		COMMON X(100),Y(100),N,X1(512),Y1(512),S(100),K
20205		DATA SZ/1002001./
20210		SS=S(K)**2
20240		SSS=S(K+1)**2
20270		IF(SS.NE.SZ.AND.SSS.NE.SZ)RETURN
20400		IF(SS.EQ.SZ.AND.X(K).LT.X(K+1).AND.
20500	     1	Y(K).GT.Y(K+1)) S(K)=-1001.
20600		IF(S(K)**2.EQ.SZ.AND.X(K).GT.X(K+1).AND.
20700	     1	Y(K+1).GT.Y(K)) S(K)=-1001.
20800		IF(SSS.EQ.SZ.AND.X(K).LT.X(K+1).AND.
20900	     1	Y(K).GT.Y(K+1)) S(K+1)=-1001.
21000		IF(S(K+1)**2.EQ.SZ.AND.X(K).GT.X(K+1).AND.
21100	     1	Y(K).LT.Y(K+1)) S(K+1)=-1001.
21200	100	RETURN
21300		END
21400	
22000		FUNCTION SLM2(L)
22100		COMMON X(100),Y(100),N,X1(512),Y1(512),S(100),K
22175		M=K+L
22187		J=M+1
22200		SLM2=(Y(J)-Y(M))/(X(J)-X(M))
22300		END
22400	
22600		FUNCTION SLM3(L)
22700		COMMON/SLM/A(2,200),B(2,200),C(2,200),I(200),I1(200),
22750		1 I2(200),KT,TT
22800		SLM3=A(L,KT)*TT**2+B(L,KT)*TT+C(L,KT)
22900		END